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Abstract 



In this paper, we develop a simulation-based approach to optimisation with multi- 
modal functions using slice sampling. Our method specifies the objective function as 
an energy potential in a Boltzmann distribution and then we use auxiliary exponen- 
tial slice variables to provide samples for a variety of energy levels. Our slice sampler 
draws uniformly over the augmented slice region. We identify the global modes by 
projecting the path of the chain back to the underlying space. Four standard test 
functions are used to illustrate the methodology: Rosenbrock, Himmelblau, Rast- 
rigin, and Shubert. These functions demonstrate the flexibility of our approach as 
they include functions with long ridges (Rosenbrock), multi-modality (Himmelblau, 
Shubert) and many local modes dominated by one global (Rastrigin). The methods 
described here are implemented in the R package McmcOpt. 
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1 Introduction 



Multi-modal objective functions pose difficulties for local search and derivative based meth- 
ods. Our simulation-based approach exploits a well-known duality between functional op- 
timisation and sampling to find the mode of the Boltzmann distribution with an energy 
potential specified by the objective function of interest. We exploit auxiliary slice variables 
to augment the state space and develop a Markov chain which samples uniformly on the 
slice region whilst traversing the desired modes in the original space. A simulation-based 
approach has the advantage of being derivative free and thus avoids some of the problems 
associated with optimisation of multi-modal functions. To illustrate our methodology, we 
test four standard global functions from the optimisation literature: the Rosenbrock, Him- 
melblau, Rastrigin, and Shubert. 

Our approach builds on the seminal work of Pincus (1968, 1970), Geman and Geman 
(1985), and Geman (1990) who proposed simulation-based Metropohs-Hastings (MH) algo- 
rithms for Boltzmann distributions in the context of constrained functional optimisation. 
We use a shce samphng approach for additive functionals (Edwards and Sokal, 1988) that 
scales to high dimensions. The Boltzmann distribution contains an energy level parameter 
which we can use to perform a sensitivity analysis. Subsequent research has shown that 
direct MCMC methods can be fraught with convergence difficulties as the associated chain 
can easily become stuck in a local mode; consequentially, careful tuning of algorithms is 
generally required. 

One of the advantages of our approach is that it does not require additional tuning. For 
the algorithm to be efficient the chain has to have good mixing properties; see, Tweedie 
and Mengersen (1994), Poison (1996) and Roberts and Rosenthal (1999), Mira and Tierney 
(2002) for theoretical convergence results on slice sampling. From a practical perspective, in 
the four examples considered here, our slice sampler has remarkably good mixing properties 
as it samples uniformly over the higher dimensional slice region. In all cases, with a 
reasonable set of energy levels, we can traverse the objective functions of interest within 
thousands of MCMC draws. 

Other popular simulation-based methods range from simulated annealing (Kirkpatrick 
et al, 1983, Geman, 1990), direct and evolutionary Metropolis MCMC (Liu, et al, 2000), 
particle swarm (Kennedy and Eberhart, 1995), multi-set sampling (Leman, Chen and 
Lavine, 1999), and stochastic guided search (Gramacy and Taddy, 2010, Gramacy and 
Lee, 2010, Taddy, Lee, Gray and Griffin, 2008). For example, Janson and Middendorf 
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(2005) illustrate particle swarm methods on the Rosenbrock and Rastigin functions. Many 
of these methods require careful tuning. Our approach relies solely on the energy level as a 
free tuning parameter and requires only a sensitivity analysis to this parameter. We follow 
the simulated tempering literature by focusing on a set of pre-determined energy values 
although our method extends easily to the case of stochastic levels as used in the Wang 
Landau algorithm 

The rest of the paper is as follows. Section 2 describes our simulation-based optimisation 
procedure. Section 3 develops slice samplers for the Rosenbrock, Himmelblau, Rastrigin, 
and Shubcrt functions. We show how to calculate the slice-set for each of the functions 
in turn and develop the associated MCMC algorithms. Finally, Section 4 concludes with 
directions for future research. 

2 Simulation-based Optimisation via Slice Sampling 

The general problem we address is to find the set of minima aigmin^^;^, f (x.) , for some 
domain X for a given objective /(x). This is the standard optimisation problem. Our 
method distinguishes itself from others in that we allow for the function /(x) to be multi- 
modal. We define the set of minima by 

^min = {x e A" : /(x) = min /(y)} . 

y 

We will find Xmin by simulation, however, we do not directly simulate /(x), rather we 
exploit a well-known duality between optimisation and finding the modes of the Boltzmann 
distribution with energy potential /(x) defined by the density 

7r«(x) = exp {-Kf{x.)} /Z^ for x e 

where — exp {— «;/(x)} dx is an appropriate normalisation constant or partition 
function. Clearly the minima of /(x) correspond to the modes of 7rK(x). One advantage 
of our simulation approach is that we do not require explicit knowledge of Z^. The only 
tuning parameter is k which is an energy level parameter of the Boltzmann distribution. 

The limiting cases of the Boltzmann distribution where k — )■ or k — > oo are of 
particular interest. They both lead to a uniform measure but on different sets. When 
K ^ 0, the limiting distribution, denoted by 7ro(x), is a uniform measure on the set X. 
When K — )■ oo, the limiting distribution, denoted by 7roo(x), is a uniform measure on the 
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set of modes, Therefore, if we can sample from the Boltzmann distribution we can 

identify the minima of the original function. Specifically, we have 

lim 7r«(x) = 7roo(x) = \A:m^n\~^Sx^^^{x.) 

where 5 denotes a Dirac measure. 

Once our problem is re-written in terms of the Boltzmann distribution we can extend 
existing methods for finding the modes. For example, Pincus (1968, 1970) proposed a 
Metropohs algorithm to simulate draws. We denote the reahsation of the Markov chain 
by X^^\X^^\ . . . , X^'^\ . . . which has equihbrium distribution, n^i^x)- Then, under mild 
Harris recurrence conditions, given any starting point X^°^ and energy level k, we have the 
limit 

Jim P (X(^) e A\X^^^ = x) = 7r^{A) 

for any Borel sets A. See Tierney (1994) and Azencott (1988) for further discussion. We 
can then use the ergodic mean ^ ^^=i -^^^^ along the chain as an estimate of the mean 
and hence, in the uni-modal case as k — > oo this will find the mode. 

There are, however, many possible Markov transition dynamics that have the appropri- 
ate equilibrium distribution. Thus the main issue becomes which Markov chain to use. We 
argue for the use of slice sampling methods. The practical insight of using augmentation 
and slice sampling is that essentially we have put some volume back into the spiky multi- 
mode regions. After the chain has converged in the higher dimensional set, we can project 
the draws back down into the dimension of interest and the chain will have no difficulty in 
traversing the modes even for lower energy values. 



2.1 Slice Sampling 

In this section we describe the developments in slice sampling and then show how to 
apply them to our optimisation problem. Suppose that we wish to sample from a possibly 
high dimensional un-normalised density /t(x). We do this by sampling uniformly from the 
region that lies under the density plot of tt. This idea is formalised by letting u be an 
auxiliary "slice- variable" and defining a joint distribution 7r(x, m) that is uniform on the 
set U = {(x, u) : < M < 7r(x)}. Therefore, p(x, m) = 1/Z on U and zero otherwise. Here 
Z — 7r(x)(ix is the appropriate normalisation constant. The marginal is the desired 
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normalised density as 



TT 



(x) = f n{x,u)du^ (1/Z) f 



ciC/ = 7r(x)/Z . 



Ju Jo 



We are then left with sampling from the uniform density on U. Neal (2003) provides a 
general slice algorithm. When it is straightforward to sample from the "slice" region defined 
by u, namely 5„ = {x : u < 7r(x)}, then a simple Gibbs sampler which iterates between 
drawing a uniform (m|x) ~ f/m(0,7r(x)) and (x|m) ~ Uni{Su) provides a Markov chain 
with the joint distribution 7r(x, u) and hence we can obtain marginal draws from 7r(x)/Z. 

The Swenden-Wang algorithm (Edwards and Sokal, 1988) extends the apphcation of 
slice sampling to product functionals. Suppose that we wish to sample from a density that 
is a product of functions: 



To do this, we introduce K auxiliary uniform slice variables {ui, . . . ,uk) with a joint 
7r(x, Ml, ... , uk) defined to be uniform on the "slice" region: 

<S = {(x, m) : Ui < 7ri(x) , 1 < i < K }. 

The marginal distribution 7rx(x) = J T^{yi^u)du — 7rj(x)/Z^. We can sample this 
distribution using the complete conditionals 

(Mj|x) ~ C/m(0, 7rj(x)) for i = 1, . . . , and (x|m) ~ t/m((Su) . 

where iS„ = {x : < 7ri(x) , I <i < K }. 

We will be interested in additive objective functions: /(x) = X^^i fii^)- Define 7ri(x) = 
exp(— «;/j(x)). Now the Boltzmann distribution is 



The exponential slice sampler extends slice sampling to additive functionals by letting 
y — {yi, . . . ,yK) be a vector of exponential variables with {yi\K) ~ Exp{k). The joint 
distribution is given by: 



K 



TT/f (x) = 7ri(a;) . . . 'Kk{x)IZk = -Kii^jZK ■ 



TTKi^) = exp (-«;/(x)) /Z« = exp -k 
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The MCMC algorithm uses the complete conditional p(x|yi, . . . , y^) and p(yi|x) where 

7r(x|y) - Urn {{/,(x) > , Vz } = u£i (x, G }• 

The auxiliary variables are sampled from (|/i|x) ~ Exp{k)1{Q, fi{^)). 

So far we have constructed a Markov chain on the augmented space (x, y) which con- 
verges in distribution to 7r(x, y): we write {^^'^\y^\ . . . , i/^'' j = (x, yi, . . . , yx) ~ 7r(x, y) 
as n — > oo. Given weak convergence, we also have for any functional that 

F{^(-\y[-\...,yP)RF{^,y,,...,yx) 

Hence we can project that joint draws back to the original space and view x*^") as traversing 
the Boltzmann funcion or equivalently /(x). This allows us to traverse the modes of interest 
and has the advantage of scalability in K. 

A related approach is to allow the energy level, k, to be random. One places pseudo-prior 
weights, denoted by and simulates the mixture Boltzmann distribution ^«;P(fi;)vrK(x). 
Another line of research is based on simulated annealing (Kirkpatrick et al, 1983, Aarts 
and Korst, 1988, Van Laarhoven and Aarts, 1987) which increases the energy level with 
the length of simulation in an appropriate schedule (Gidas, 1985). Other approaches that 
randomize k include multi-canonical sampling (Berg and Neuhaus, 1982), simulated tem- 
pering (Marinari and Parisi, 1992, Geyer and Thompson, 1995) uses a random walk on a set 
of energy level, equi-energy sampling (Kou, Zhou and Wong, 2006), evolutionary MCMC 
(Liu, Liang and Wong, 2000, 2001) and the Wang-Landau algorithm (Liang, 2005, Atchade 
and Liu, 2010) and When X is bounded and supp{f) < oo, shce sampler has the added 
property of geometric convergence (Roberts and Poison, 1994) and uniformity (Mira and 
Tierney, 2002). 

3 Four Examples 

Figures 1 and 2 show contour and drape-plots of the Rosenbrock, Himmelblau, Rastigrin 
and Shubert functions. This set of functions exhibits a variety of challenges: a global mode 
in a long valley (Rosenbrock), multiple local and one global mode (Rastrigin) to multi- 
modality (Himmelblau, Shubert). For example, traditional derivative-based methods have 
difficulties for the Rosenbrock function to traverse its long steep valley. 
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Let /(x) = Z]i=i 2:2) be a bivariate additive objective function defined over a 
bounded region for x = {xi, X2} G 3?^. We will apply the exponential slice sampler as the 
functions to the Boltzmann distribution 



First, as in simulated tempering, we have to define a set of temperatures {ki < . . . < Km} to 
run our Markov chain. For all four examples, we pick m = 4, and we use /t G {0.1, 0.5, 1, 5} 
except for the Rosenbrock function where we set k e (1, 5, 50, 5000} which requires higher 
energy levels. 

3.1 Rosenbrock function 

Rosenbrock's valley is a classic optimization problem that illustrates the difficulties with 
local search methods. The global minimum lays inside a long, narrow flat valley. Finding 
the valley is straightforward; however, getting to the minimum is hard. We need to find 
the minimum (xi, X2) = (1,1) of the function: 

/(xi, X2) = (1 — xif + c{x2 — x\Y where c = 100 . 

The Boltzmann distribution has density 

t^k{xi,X2) = exp (-K {(1 - xif + c{x2 - xlf}) /Z^ . 

There are a number of ways of introdTicing slice variables. We choose to slice out the last 
nonlinear factor. Let (M|a;i,a;2) ~ Uni {0,exp {—kc{x2 — xl)^}). Then we have a three 
variable joint distribution: 

7r^{x, y,u) = exp {-k(1 - Xi)^} I (O < u < exp [-kc{x2 - xj)'^]) /Z^. 

We can implement MCMC using a partially collapsed Gibbs sampler (van Dyk and Park, 
2008, Park and van Dyk, 2009). That is, we can marginalise u out of the draw for X2 and 
use the conditional 7r(x2|xi) rather than Tr{x2\xi,u). The complete conditionals are then: 




'k{x2\xi) M{x\,2/ kc), 
Ti{xi\x2, u) ~ A/'(l, 2/ k)1 {a{u, X2) < x < b{u, X2)) , 
Ti{u\xi, X2) ^ Uni (0, exp \^—kc{x2 — a^i)^}) ■ 
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The interval {a{u, X2), b{u, X2)) is found by inverting the shce region: 

u < exp [—Kc{x2 — xl)^} imphes X2 — \J — hiu/ kc < x\ < X2 -\- y/— Ihu/kc. 
For 6 > and a < xf < b, we have —^/a < xi < Vb and so 

a{u, X2) — —\Jx2 — -x/— Iuu/kc and b{u, X2) — \J X2-\- \l — Iuu/kc. 

Figure 3 shows a sensitivity analysis for a range of energy values k G {1,5, 50, 5000}. We 
run our MCMC algorithm for G = 1000 with a burn-in of Go = 100. Higher energy levels 
are required for the chain to travel along the valley to the minimum at {xi,X2) = (1,1)- 
As we increase k, the shce sampler is able to traverse the valley and find the minimum. 

3.2 Himmelblau's function 

Himmelblau's function is defined by 

f{xi, X2) = {xl + X2- liy + (xi +xl- 5)1 

This function has four identical local minima at zero and a local maximum at {xi,X2) — 
(-0.27, -0.92). The minima are at 

(3, 2), (-2.805, 3.131), (-3.779, -3.282), (3.584, -1.848) 

with a function value of zero. The associated Boltzmann distribution is 

t^k{xi-,X2) = exp{-«; [{xl + X2 - llf + {xi + xl - bf)] /Z^ . 

Due to the quadratic terms also containing squares this distribution is not straightforward 
to simulate from. We observe that the following inequalities hold for the minima: x\-\-X2 — 
11 < and xi + - 5 > 0. 

To implement slice sampling, we use a latent variable augmentation u — (wi, ^2) and a 
joint distribution 

tIk{xi, X2, Ui, U2) = I (0 < Ml < exp {^—i<i{x\ + X2 — 11)^}) I (0 < M2 < exp \^—k{xi -\- x\ — 5)^}) /Z,^ . 
Given (1^1,^2), we can invert the slice regions to obtain the inequalities: 

—Kr^\ogui>{x\-\-X2 — ^^Y and — \0gu2 > {xi -\- x\ — . 
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Therefore, for {xi\x2), we have 

ai = 11 — X2 — a/— log Ml < x1 < 11 — X2 + sj —K^^ log Ux = hi] 



Cl 



5 — xl — \J logM2 <X\<h — x\ + logM2 = (il- 



Given the inequahties: < x\ < hi and Ci < Xi < di, we can first without loss of 
generality replace Oi — > max(ai,0) and assume that Oi > 0. Then, we have the union of 
the following regions: 



— V 6i < xi < —^/oi and y/oi < xi < hi. 

Combining, we have 

max ^— a/^i, ci^ < xi < min (— -^/oi, di) or max (-v/oi, ci) < xi < min ^-xAi; ^i) 

For the conditional 7r(a;2|a;i) we can argue in a similar fashion to obtain the constraints 

02 = 5 — xi — \/ —K~^ log < X2 < 5 — xi + y^—K~^ log = 62; 
C2 = 11 - - -xZ-K-Uog ^2 < a:2 < 11 - Xi + ■\/—K-'^ log ^2 = (^2- 

The complete set of conditionals is then given by: 

7r(a;i|x2, Ml, M2) ~ t^?^^ ^max ^— -^|ai|, Ci^ ,min ^\/&i, rfij j , 

7r(x2|xi,iti,M2) ~ C/m ^max ^-•\/|«2l, C2^ ,min [\/h2,d-^^ , 

7r(Mi|a;i, X2) ~ Uni (O, exp |— /t(xi +X2 — H)^}) , 
7r(M2|a;i, a;2) ^ Uni (O, exp |— «;(a;i + — 5)^}) . 

Figure 4 shows a sensitivity analysis to k, G {0.1, 0.5, 1,5}. The slice sampler is again run 
for only G — 1000 iterations with a burn-in of Go — 100. With longer chains and larger 
energy levels the algorithm will traverse the four modes with equal probability. Of the 
examples that we consider here, the Himmleblau function would benefit the most from a 
mixture energy level distribution to traverse the contours of the underlying function. 

3.3 Rastrigin 

The 2-dimensional Rastrigin function is defined on the region —5.12 < xj < 5.12 by: 

2 

f{xi, X2) = 2A + ^ (x| - A cos(27ra;j)) with A = 10 
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with a global minimum at {xi, X2) — (0, 0). The Boltzmann distribution then becomes 

T^Kixi, X2) = exp I -K ^ ,1;^ I cxp I kA ^ cos(27rXj) | /Z^ . 

We use exponential slice variables {yi, 1/2) and a joint distribution defined by 

7r«(xi,X2,yi,y2) = exp |-/t^a;|| I(-K^cos(27rXj) < yj) e'^'/Z^ . 

The slice region is invertible and is specified by the set of inequalities for j = 1, 2 

cos(27rXj) > {—yj/AK,) . 

The conditional Xj draw then results from a normal draw restricted to this interval set. A 
Gibbs sampler with the conditionals for j = 1, 2 is as follows: 

j, y) r^Af (0, {2k)-^) I (cos(27rXj) > {-yj/AK)) 
TrK{yj\xj) ~ Exp[— K A cos{27rxj), 00) 

For the implementation over Xj G [—5.12, 5.12], we draw the truncated normals and trun- 
cated exponentials for the shce variables. 

Figure 5 shows a sensitivity analysis to k e {0.1, 0.5, 1, 5} with G — 1000 and a burn-in 
of Go = 100. Again slice sampling of the Boltzmann distribution finds the mode in a 
straightforward manner. 

3.4 Shubert function 

The Shubert function is defined within the domain I(— 10, 10) by 

5 

f{xi,X2) = —C{xi)C{x2) where C(a;) = j cos ((j + l)a; + j) . 
We need to simulate from the Boltzmann distribution 

For the conditional TrK{xi\x2), we can write 

5 

7^k(xi\x2) = J]e"^(^^)^'^°^((^'+^)^i+^'). 
10 



We introduce a set of auxiliary slice variables yj,l < j < 5 for each Xi,i — 1,2 that are 
conditional exponentials. The corresponding joint is: 

5 

7r«(xi, X2, yi, ■ ■ ■ , ys) = ex.p{-K,yj)I {yj > C{x2)j cos {{j + l)xi + j))) /Z^ 
This is inverted using the condition: 

This gives a collection of intervals Ii{yj,X2) for 1 < j < 5 for each of the shce variables; Xi 
is then uniformly distributed over the intersection of these intervals, r\jIi{uj,X2). 
Hence we can then run a Gibbs sampler with the conditionals, for 1 < j < 5: 

7iK{xi\x2,y) ~ Uni{n^j^Ji{yj,X2)), 

7r^{yj\xi) Exp [C{x2)j cos {{j + l)xi + j) , oo) . 

Similarly, this process is repeated for 7ri^{x2\xi). This defines a 12-dimensional Gibbs sam- 
pler that is able to traverse the joint distribution. 

Figure 6 shows a sensitivity analysis to the same set of energy levels used for the 
Rastgrin and Himmleblau functions, namely k G {0.1, 0.5, 1, 5}. Again the projected draws 
from the uniform slice region traverse the modes of the associated Boltzmann distribution 
in an efficient manner. 

4 Discussion 

We have described how slice sampling methods can be applied to functional optimisa- 
tion. Our approach is parallelisable as in shce sampling (Tibbits et al, 2011). While 
we have only considered four test functions our methodology applies to a multitude of 
multi-modal functions, see Molga and Smutnicki (2005) for a hst of possible candidates. 
Our approach is flexible enough to handle additive functions that are multiplicative. For 
example, the Michalewicz function is defined by f[xi,X2) = — ^^^^ sm(a;j)sm^™ (ixf/vr) 
with m = 10 (see Yang, 2010a,b). This function has a minimum of —1.801 at the point 
{xi,X2) = (2.20319,1.57049) and it is challenge to fix the minimum. We also note that 
certain functions are straightforward as the Boltzmann distribution is conditionally normal. 
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For example, the Booth function defined by /(xi, X2) = {xi + 2x2 — 7)^ + {2xi + X2 — 5)^ has 
a minimum of zero at the point {xi,X2) — (1, 3). The minimum can be identified without 
resorting to simulation as the Booth function can be factorised into the quadratic form 
(x — ii)'Q~^{x — jj) where x — (xi, X2) with jjt, — (1, 3) and Q — \ {^, 4; —4, 5). Therefore, 
the minimum can be immediately identified without resorting to simulation. 

There are clearly many other applications of these methods. For example, simulated 
annealing methods have been proposed in mixed integer non-linear programming problems 
(Cardoso et al, 1997) and in constrained optimisation (Gcman and Geman, 1985, Geman, 
1990, Whittle, 1992, Birge and Louveaux, 1997, Mueller, 2000, Asmussen and Glynn, 2008). 
Slice sampling the Boltzmann distribution provides a flexible alternative. 
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Figure 2: Drape plots: Himmelblau, Rastigrin, Shubert and Rosenbrock 
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Rosenbrock: kappa=1 



Rosenbrock: kappa=5 
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Figure 3: Rosenbrock function: f{xi,X2) = (1 — XiY + 100(a;2 — xlY- The minimum occurs 
at {xi,X2) = (1, 1). 
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Himmelblau: kappa=0.1 



Himmelblau: kappa=0.5 




Figure 4: Himmelblau function: f{xi,X2) = {xf + X2 — 11)^ + (xi + X2 — 5)^. There are 
four idenital local minima at zero and a local maximum at Xi = — 0.27, X2 = —0.92. 
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Rastigrin: kappa=0.1 



Rastigrin: kappa=0.5 
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Figure 5: 2-d Rastigrin function: f{xi,X2) = 2A + Y^'^j=i{x'j — Acos{2TTXj) with A = 10 and 



-5.12 < Xj < 5.12. 
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Shubert: kappa=0.1 
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Figure 


6: Shubert function: f{xi,X2) = 


-C(X1)C(X2) 


where 


C{x) 



Yl]=i3cos{{j + l)x + j). 
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